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A calculation of the interaction potential of two heavy-light mesons in lattice QCD is 
used to study the existence of tetraquark bound states. The interaction potential of the 
tetraquark system is calculated on the lattice with 2+1 flavours of dynamical fermions with 
lattice interpolating fields constructed using colorwave propagators. These propagators pro- 
vide a new method for constructing all-to-all spatially smeared the interpolating fields, a 
technique which allows for a better overlap with the ground state wavefunction as well as 
reduced statistical noise. Potentials are extracted for 24 distinct channels, and are fit with a 
phenomenological non-relativistic quark model potential, from which a determination of the 
existence of bound states is made via numerical solution of the two body radial Schrodinger 
equation. 



I. INTRODUCTION 

The calculation of hadronic forces from first principles allows insight into how interactions of 
the fundamental quark and gluonic degrees of freedom manifest themselves at the hadronic level. 
Lattice QCD is an excellent tool for calculating hadronic observables in the low energy regime. 
Although lattice calculations in euclidean space are not well suited for the study of real-time 
scattering processes, two methods can be used to extract interaction information from the lattice. 
One method, developed by Liischer [1], relates the elastic scattering phase shift of a two particle 
system in a finite periodic box with the energy levels of the system. An alternate method, used in 
the present work, extracts the interaction energy as a function of hadron separation. This method 
is only applicable for systems of hadrons containing more than one heavy quarks which can be 
treated in the static approximation providing a definite spatial position for the hadrons. 

Phenomenologically, two heavy-light meson systems (which we will denote as HLHL) have 
become interesting in the study of tetraquark bound states [2] [2] [I]- It has long been known 
that the binding of a QQqq (with q = u,d) system increases with the mass ratio of the heavy 
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to light quark flavours [5], thus ccqq and hhqq systems are excellent candidates in the search 
for exotic four quark bound states. In Ref. [1] a distinction was made between two types of 
tetraquark bound states: molecular, in which the four quarks exhibit a single physical two- meson 
(singlet-singlet) component, and the more exotic compact bound states. The latter would involve 
a complicated color space structure in which quark pairs form color vectors which then combine to 
form a colorless four quark state [4J. In spite of this complicated color structure, compact bound 
states can be interpreted as a mixture of various two meson (color singlet) components [6j. The 
expected features that would characterize a molecular bound state would be a small binding energy 
and a bound state RMS radius greater than that of the sum of the two particle sizes, i.e.: 

RMS,, 
^ RMSm^ + RMSm2 

A compact state, on the other hand, would be more tightly bound and have a smaller RMS 
radius than the molecular state. In Ref. [7] doubly heavy four quark states were modeled as 
hadronic molecules interacting via a meson exchange potential. Several of the doubly bottom bound 
states were found to be deeply bound and spatially compact, making them excellent candidates 
for tetraquark bound states. It is with these ideas in mind that we may begin to search for the 
signature of compact bound states on the lattice. 

A recent lattice calculation of the HLHL interaction energy [8j has in fact hinted at the possibility 
of a bound tetraquark state in one channel that exhibits a particularly wide and deep potential 
well when compared with other channels, although no exhaustive determination of a bound state 
was undertaken. Our work goes beyond this presenting a quantitative determination of a bound 
state energy in the HLHL system from a lattice calculation. 

An inherent difficulty in making comparisons between theoretical models and lattice calculations 
performed in the static limit stems from the omission of the heavy quark spin in the static limit. 
As rriH — >■ oo, the integer valued ( J = 0, 1) angular momemtum eigenstates of a single heavy light 
meson map onto a single static limit eigenstate with J = 1/2. The energies of the non-static angular 
momentum eigenstates also converge to a single energy corresponding to the J = 1/2 eigenstate. 
Although the two spaces map onto each other, there is not a simple one to one correspondence 
between static limit eigenstates and their non-static counterparts, and care must be taken in making 
identifications between the two spaces. Previous lattice studies of HLHL interaction energy ([S], 
|10j for example) performed in the quenched approximation and included uncontrolled systematic 
errors because of this. Recently dynamical quarks have been used to calculate the HLHL interaction 
energy using a complete set of quantum numbers which exploits the full set of symmetries of the 
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HLHL system p]. 

With our choice of quantum numbers (presented in section |ll]) we are able to draw a connection 
between the quantum numbers and the quahtative behavior of the states. Additionally, by way of 
symmetry arguments, we are able to relate our static- limit states to non-static angular momentum 
eigenstates. 

II. BACKGROUND 
A. Heavy-Light states 

The quark model view of a heavy-light meson is of a heavy anti-quark Q coupled to a light 
quark q. The relevant quantum numbers to describe such a state are total angular momentum J 
and its projection along some axis (here arbitrarily chosen to be z) J^, and the parity Pi as well as 
the relevant flavor quantum numbers. For our interests, we choose Q = b and q = {u, d}. Therefore 
all states then have bottomness 6 = +1, and are otherwise classified by total isospin and the third 
component of isospin {I,Iz) = (1/2, ±1/2). Throughout this work, we make the assumption that 
we fit our correlation functions with a sufficiently large tmin such that contributions from excited 
states have died out and we extract only the ground state energy. Furthermore, we assume that 
states with non-zero orbital angular momentum L are at sufficiently high energies as to have a 
negligible contribution to the ground state energies which we extract. We are then free to speak 
of the spin and total angular momentum interchangeably. 

In heavy quark effective theory, spin dependent contributions enter into the heavy quark action 
at order l/niH, and in the static limit {mn oo) the heavy quark acts as a static color source. 
This means that the spin of the HL meson comes only from the light degrees of freedom. Because 
of this, the physical HL meson states with J = (0, 1) become degenerate in the static limit, with 
both represented by a single J = 1/2 state. The relevant angular momentum classification is then 
{J, Jz) = (1/2, ±1/2). With the above flavor assignments, the lowest energy excitations of the 
B spectrum with = {0, 1}~ (coupling to the static = l/2~ B) are -Bo,± and B* , and for 
= {0,1}"*" (coupling to the static = 1/2+ Bi), the ground state i?i (5721)° (neglecting 
excited states). 
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B. Heavy-Light Heavy-Light states 

When constructing states with a pair of HL mesons, care must be taken in determining a relevant 
set of quantum numbers that fully exploit the symmetries of the problem. The flavor quantum 
numbers for a Heavy-Light Heavy-Light (HLHL) system are straightforward, and for a QqQq (with 
q = {u, d}) there are two isospin combinations, an isospin triplet with 7 = 1 and an 7 = singlet. 
For a HLHL pair separated by a vector r the rotational symmetry is broken to rotations around the 
separation axis. Total angular momentum J is therefore no longer a conserved quantity, though its 
projection along the axis of separation (arbitrarily taken to be z) is still conserved. The system will 
also be symmetric or antisymmetric under parity as well as reflections through a plane containing 
the separation axis, which we shall call Pj^. This last transformation can be accomplished by a 
parity transformation followed by a rotation of tt about an axis perpendicular to the reflection 
plane. States with J2 = ±1 are not invariant under this transformation (being mapped onto 
each other), but their average is an eigenstate of P±. Lastly we choose to classify HLHL states 
by intrinsic parity Pi, defined to be the product of the intrinsic parities of the two light quarks, 
and (full) parity P, defined as the product of the intrinsic parity transformation and coordinate 
inversion of the two particle spatial wavefunction. We will use both parity quantum numbers in 
our classification of states. 

III. METHODOLOGY 
A. HL and HLHL interpolating fields 

A general interpolating operator coupling to a single heavy-light state is given by: 

OHL{x)=Qix)rqix) (1) 

with r chosen to achieve the desired angular momentum and parity quantum numbers. For pseu- 
doscalar HL states, F = 75, 7j (with z = 1, 2, 3), corresponding to a particle in the static limit with 
= l/2~, which we will refer to simply as B. J = 1 meson states with F = 1,7^75 correspond 
to a state with = 1/2"'", which we shall refer to as Bi. We make the choice F = 75 for Ob and 
T = 1 for Obi- As it will be useful in the analysis of HLHL states, it should be noted that for these 
choices of F, correlation functions constructed from Ob interpolating fields will consist of only 
upper (positive parity) components in the Dirac basis of the light quarks while those constructed 
from Obi will consist of only lower (negative parity) components. This is explicitly shown in Ap- 
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pendix [A| The states are classified by the additional flavor quantum numbers (I, Iz) = (1, ±1) for 
q = {u,d}. 

For HLHL states, we want to create states with definite (/, Iz, \Jz\, P±, P, Pi) and displacement 
r at the source and sink. To do this, we want to couple only our light quarks in spinor space to 
specify the quantum numbers of the state while allowing the heavy quarks to act only as color 
sources. Our general HLHL operator is then given by: 



^SJl^''''^'''''''^ ^ = Q t)Q{x + f, t) X [q (x, t)q{x + f, t)] 



(2) 

il,h,\j,\,P^,P,Pi) 



where the light quark wavefunctions [q {x, t) q{x + f, t)] are combined in such a way as to achieve 
the set of quantum numbers (/, Iz,\Jz\i P±, P, Pi) of the system. The explicit construction of these 
wavefunctions is described in Appendix |B] For simplicity we restrict ourselves to identical source 
and sink interpolating fields neglecting any cross correlators between states. Isospin is a good 
quantum number on the 2 + 1 flavor lattices with which we work, and we choose our interpolating 
fields to be isospin eigenstates with (/, Iz) = (1, 1) and (/, Iz) = (0, 0). At large spatial separations, 
we expect the energy of the four quark state to asymptotically approach the energy of it's dominant 
two meson component^. States with Pj = — 1 will tend towards the energy of a BBi combination at 
large spatial separations. There are two possible combinations of the light quark parities that yield 
Pi = +1: {pi,P2) = (+, +) , {—, —)■ In light of the above discussion of parity content of single HL 
states, we project our Pj = +1 interpolating fields to contain only negative or positive parity spinor 
components and retain these as distinct interpolating fields. The expectation is that interpolating 
fields constructed from lower spinor components will exhibit a significantly higher ground state 
energy in relation to those constructed from upper components. The reason for this is that the 
(— , — ) interpolating field are constructed by the product of two Pi meson interpolating fields, thus 
should exhibit an asymptotic energy (as r — t- oo) near twice that of the single Pi energy. Similarly 
the (+,+) interpolating field is constructed from the product of two P meson interpolating fields 
tending asymptotically as f — )• oo towards a ground state energy of twice that of a single P meson. 
We differentiate all interpolating fields by their dominant asymptotic content in the tabulation of 
interpolating fields in Table |I| 



^ Here we are referring to the dominant lowest energy contribution, as we expect excited states to contribute 
negligibly to the extracted HLHL ground state energies 



6 







{1 7 IzAJzl tP-L 


P,P^) 


Dominant asymptotic content 


(1,1,1,-,- 


-,+) 


(0,0,1,-,+,+) 




BB 


(1,1,0,-,- 


-,+) 


(0,0,0,-,+,+) 




BB 


(1,1,0,+, 


+,+) 


(0,0,0,+,-,+) 




BB 


(1,1,1,-,- 


-,+) 


(0,0,1,-,+,+) 




BiB, 


(1,1,0,-,- 


-,+) 


(0,0,0,-,+,+) 




BiBi 


(1,1,0,+, 


+,+) 


(0,0,0,+,-,+) 




BiBi 


(1,1,1,+, 


+,-) 


(0,0,1,+,-,-) 




BBi 


(1,1,0,+, 


+,-) 


(0,0,0,+,-,-) 




BBi 


(1,1,0,-,- 


-,-) 


(0,0,0,-,+,-) 




BBi 


(1,1,1,+, 


-,-) 


(0,0,1,+,+,-) 




BBi 


(1,1,0,+, 


-,-) 


(0,0,0,+,+,-) 




BBi 


(1,1,0,-,- 


K-) 


(0,0,0,-,-,-) 




BBi 



TABLE I: HLHL interpolating operator basis and expected asymptotic values 



IV. DETAILS OF THE LATTICE CALCULATION 

We work with colorwave propagators (described below) calculated on = 2 + 1 anisotropic 
(24^ X 128) lattices generated by the Hadron Spectrum Collaboration |12j with a pion mass of 
roughly 380 MeV. The fermion action used was the clover Wilson action with stout link smearing, 
not smeared in the temporal direction. The gauge action was Symanzik tree level tadpole-improved 
without a rectangle in the temporal direction, preserving temporal ultra-locality. The spatial and 
temporal lattice spacings for these lattices are as = 0.1227(8)fm and at = 0.03506(23)fm. The pion 
mass on this ensemble is 0.0681(4) in temporal lattice units. The Chroma Software package for 
Lattice QCD [13J was used to generate both colorwave and heavy propagators. The calculation 
of the HL and HLHL energies was performed using 305 gauge field configurations with eight 
sources spaced evenly in the temporal direction. Ground state energies were extracted using single 
exponential correlated fits, with an appropriate tmin determined from the quality of the fit. 

A. Colorwave Formalism 

1. Two quark states 



Consider a general operator for a two quark mesonic state: 

O (f ) = ^1 (f ) Tq2 (f ) 



(3) 
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where we assume for simplicity that the two quarks have different flavors. We seek to calculate 
the correlation function with localized interpolating fields: (averaged over spatial source and sink 
locations to increase statistics) 



C{t,to) = ^(0(2/)0t(i 

x,y 

= ^^tr{Si{x,to\y,t)rS2{y,t\x,to)r) (4) 

X y 

Following the methodology presented in ^4j, we now consider any complete set of orthonormal 
states {(pi (x)} which satisfy: 

^ 0* (x) {y) = 6{x-y),Y^ 4>* (x) 4>j (x) = S,j . (5) 

i X 

By inserting the completeness relation of eq. [5] twice into the two point function of eq. |4j 
C (i, to)=T.Yl (-^i *o|y'' t)5{y- y') TS2 (y, t\x',o ) S {x - x') T) 

x,x' y,y' 

= E E (^1 E (y) (y') {y, t\x', to) <p* (x) (x') r\ 

x,x' y,y' \ i j I 

= Y^Srito,t)TS'^^it,to)T (6) 
where we have defined: 

5^'^' it, to) (y) ^ <^j- 

x,y 

A convenient choice for the {4>i{x)} is a plane wave basis: (j)i{x) = 4>p{x) = e~'^'^^5s^s'^c,c' ■ 
The delta functions here operate on color and spin. With this choice of basis, we define 
S^'^t,to) = SP'P {t,to) to be colorwave propagators. The use of these propagators allows us 
to implement spatial smearing at the source and sink of our correlation functions. In the limit 
where all momenta are summed over in equation [Tj all to all point-point propagators are recovered. 
However, introducing a maximum momentum cutoff p^^^ we are able to introduce and control the 
effective amount of spatial smearing^. The effect of restricting the plane wave basis to |pp < p^^j 
(summing over a momentum space volume) is illustrated in Fig. [l] where effective masses for 
single HL B meson correlation functions^ are presented. It's evident that the noise of the signal 



^ It should be noted that smearing is achieved only in fixed gauge. In our case we use the Coulomb gauge, which is 
a smooth gauge allowing to project out high energy modes if the cutoff p^^t is kept small . 



These IfL correlation functions are defined in Appendix |Aj eq. Af 
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decreases by increasing the momentum space cutoff (as tliis increases tlie statistics contributing to 
the correlation function). 

Each effective mass plateau appears to begin at roughly the same point independent of p^ut^ 
and thus a common fit range of 17 — 30 was chosen for all values of p^^^. In Fig. [l]we can see that 
as p^^f increases the overlap with excited states drops resulting lower values for the effective mass 
at earlier times. This indicates that a small radial smearing of the quarks field results interpolating 
fields that have better overlap with the ground state of the system. Such behavior is likely due to 
the fact that the a non-relativistic HL meson in the static limit is a highly localized object whose 
wavefunction is confined to a small spatial region. 

In light of this behavior and in order to reduce computational cost associated with increasing 
the momentum cutoff, a value of = 1 was chosen for calculations of the HLHL system. 

2. HLHL States 

We begin with a correlation function for two heavy-light mesons separated by f as described 
above: 



Each heavy quark source can only be contracted with the sink at the same spatial location, 
and upon contraction we work only with the Wilson line portion of the heavy quark propagator, 
as we want the quantum numbers of the system to be determined entirely by the light degrees of 
freedom. There are two possible light quark contractions, one where the light quarks contract with 
source and sink at the same spatial location (direct), and one where the light quarks contract at 
the other spatial location (crossed). Performing these contractions, we have (omitting the overall 
color trace): 




X 



= ^ {Q {x, t)Q{x + f, t) q (x, t)q{x + f, t)q{x + f, to) q (f , to) Q {x + f, to) Q (x, to)) 



X 



Chlhl (t, r) =Y^ 75 W^"^ (x; t, to) 7575^"^ (x + f; t, to) 75 



X 



X tvd [S {x + f,t;x + f, to) "S" (x, t; x, to) - 5 (x + f, t; x, to) S (x, t; x + f, to)] (9) 



Here, tr^ denotes the trace over Dirac space spinor indices and W is the Wilson line 



t 




(10) 
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FIG. 1: Effective mass for HL B for increasing \p^\ < \pf 
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We now introduce our partially fourier transformed light quark propagators as: 

S{x[,t;xi,to) = ^ e^P'i<S{p[,t;p,,to)e-'P^-^ (11) 
p'l ,Pi 

where sums over momenta pi have been restricted to < 1 as described in the previous section. 
Using this, the above correlator can be rewritten as: 

Chlhl {t, ^= H ^^^^ ^575 {x + r; t, to) 75 x e^iPi-n+p'^-P^^e^ip'^-P^y 

X [S {p2,t;p2,to) S {pi,t;pi,to) - S {p2,t;pi,to) S (p'l, t;p2, to)] (12) 

Defining 

V (f, t, to, u;) = Y, 75 (x; t, to) 7575 {x + f; to, t) 75e^(")" (13) 
with = p'l — Pi + P2 — P2, our the final form of our HLHL correlation function becomes: 

PlP'lP2P'2 

X [S (p2,t;p2,to) S {p'i,t;pi,to) - S {p2,t;pi,to) S {p'i,t;p2,to)] (14) 

With this method, we calculate the costly D {f,t,tQ,uj) first using a parallel code (paralleliza- 
tion over space time) and then perform the far less expensive contractions with the colorwave 
propagators for our complete operator basis on a scalar workstation class machine. 

V. HLHL RESULTS 

For q = {u, d} we have 24 unique HLHL corresponding to the operators enumerated in Table 
m Each potential curve is calculated by taking the jackknife difference between the energy of the 
HLHL state for various f and the energy of the expected two meson asymptotic state: 

V (f) = Ehlhl (r) - Eb,,, - Eb,,, (15) 

The statistical uncertainty for each point is determined from jackknife statistical analysis. The 
systematic uncertainties are determined by adjusting the chosen fit range by one time slice in each 
direction and averaging the observed deviations in the energy. The systematic uncertainty for both 
Ehlhl and EB^^-^^ are determined independently and then added in quadrature to determine the 
systematic uncertainty on V (r). 



11 



We find three different asymptotic values for the various states as illustrated in Fig. [2]. The 
lowest lying asymptotic value corresponds to states with a positive intrinsic parity Pi with all spin 
components in the correlation function projected to the upper spin components, while the highest 
asymptotic value corresponds to states with positive intrinsic parity and all spins projected to the 
lower components. This asymptotic behavior is in line with our expectation that the spin projection 
of our positive intrinsic parity operators helps to increase the coupling to the lower energy BB state 
or the higher energy BiBi state. The energy difference between the highest and lowest asymptotic 
values is roughly twice the energy difference between the single HL B and Bi states, indicating 
that they are both tending asymptotically towards their expected two meson asymptotic energies 
at long distances. The slight overshoot of the highest asymptotic state beyond it's expected value 
of twice the Bi energy for d > 0.8 fm may be indicative of contamination from mixing of the HL 
Bi with a IT — B state. All Pi = (— ) states exhibit an asymptotic tendency towards the sum of the 
single HL B and Bi energies as expected. 

As the states with the lowest asymptotic energy values trend most cleanly towards their expected 
asymptotic value (indicating the least contamination from excited states), we will focus mainly on 
these states which we present in Fig. [s]. 

Several aspects of these potential curves should be noted: First, we find that the product 
of exchange parity P and intrinsic parity Pi, which is the symmetry of the two meson spatial 
wavefunction under spatial inversion, directly corresponds to the attractiveness (— ) or repulsiveness 

(+) of the state. This is in agreement with j8j. Second, the (/, I^, \ Jz\, P±, P, Pi) = (0,0,0) H h 

exhibits a significantly deeper and wider potential well when compared with the two other attractive 
channels. This qualitative difference was acknowledged in [8], and the quantum numbers of this 
channel are consistent with a bound state predicted in a phenomenological model in |4|. 

VI. BOUND STATES 

As the HLHL system has been predicted to be an excellent candidate for bound tetraquark 
states, we seek a quantitative method for extracting such a bound state (if one exists) from our lat- 
tice calculation. Our method is as follows: We fit our lattice potential to a phenomenological quark 
model potential as described in [15j. We make the choice to focus on the {I,Iz, \Jz\, P±, P, Pi) = 
(0, 0, 0, +, — , +) channel, as previous work has hinted at the possibility of a bound state here. As a 
control, we also perform the fit for the (/, Iz, \ Jz\, P±, P, Pi) = (1, 1, 0, — , — , +) channel as well. In 
our fit, we neglect the r = points as the finite value of the potential at r = is a lattice artifact 
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FIG. 2: Comparison of BB vs. BBi (left) and BB vs BiBi (right) asymptotic states. Here we take the 
energy difference for the three potential curves with respect to twice the HL B energy 



stemming from the ultraviolet cut off introduced by the lattice discretization, leaving us with 7 
data points for each potential curve, and two free parameters from the fit model. The model with 
the extracted fit parameters is then taken to be the interaction potential between two B mesons 
in the continuum limit. The two body (one-dimensional) Schrodinger equation is then solved nu- 
merically with this interaction potential to determine the existence of any negative energy (bound) 
states. It should be noted here that the solutions to the Schrodinger equation will converge to 
their continuum values as the continuum limit of the lattice calculation is taken. As we have only 
a single lattice spacing available to work with this continuum extrapolation is not an option, and 
it should be understood that the results presented in this section are at finite lattice spacing. 



A. Potential Model 



We have limited our displacements |r| < 1.27 fm, therefore long range effective interactions due 
to meson exchange do not provide a good description of the HLHL system. In reference [15j . a 
quark model picture of a two meson interaction was used to derive an interaction potential for the 
HLHL system, which included color coulomb, spin-spin, linear confinement interactions. Details 
of the derivation of the potential model can be found in the aforementioned reference, and we will 
only highlight several modifications we make when fitting this potential model to our numerical 
results. The quark model HLHL potential has the form: 



Vbbds (r) = CiVcc {as, 13, r) + C^.^Vss (os, /?, m, r) + CiVi^ {b, /3, r) 



(16) 
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FIG. 3: Calculated HLHL BB energies with expected asymptotic value (twice the calculated HL B mass) 
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with: 



Vcciois,f3,r) 
Vss {as,P,fh, r) 



9r 



1/2 



1/2 



I3r - 4Erf 



(?) 



-/32r2/2 



27 Vtt/ fn? 



(17) 
(18) 



3^ 



^,e-^V/2^2(2^j 



1/2 



/3r + — 1 Erf 
pr 



/32rV2 



1/2^ 



-3/32r2/4 



(19) 



Here, is the strong coupling constant, /3 is the spatial width of the quark model single HL 



meson wavefunction, m is the mass of the light quark in the MS scheme, and h is the QCD string 
tension. The coefficients Cj and Cs-s, which contain the spin information of the HLHL state, are 
defined as matrix elements between initial (unprimed) and final (primed) two meson states and will 
be discussed further below. It should be noted that the above potential model acquires an overall 
minus sign if the isospin wavefunction of the two meson state is antisymmetric. Additionally, the 
potential is a function of |r| and not r, as any tensor interaction terms are neglected in this model. 



B. Fit Model 



When applying the above model to our lattice data, we must make several modifications to the 
above quark model potential. Due to the use of periodic boundary conditions in the calculation, 
interactions with image "charges" lying past the boundary must be accounted for. We must 
also consider the possibility that there will be long range meson exchange interactions that were 

neglected in our choice of potential model. To account for these long range interactions, we extend 
the original model by adding a simple Yukawa like term for one pion exchange: 

p-m-„r 

V^-^ (r) = Vbbds (r) + (20) 

Here we take to be the mass of the pion on the gauge field configurations used in the calculation 
(~ 390 GeV). The parameter g is discussed below. 

In principle, interactions with each of the infinitely many image charges contribute to the 

potential and must be included. In practice however, we may restrict ourselves to contributions 
where the image of the first meson is < 3L/2 (~ 4.5 fm) away from the second and vice versa. This 
approximation is valid as the contribution of these truncated images (at separations of r > 3L/2) to 
the potential (with the choice of parameters outlined below) is O (lO~^) MeV. With the inclusion 
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FIG. 4: Contribution of image charges to the potential (left) and contributions to the potential model 
Vhlhl from the individual terms in eq. [16] 



of the image charges our potential model then becomes: 



^Im — ^ 



Yuk 



(21) 

ri<3L/2 

The addition of these image charges modify the potential at long distance as illustrated in Fig. [4] 
The final modification made to the potential model is a modification of the spin dependent 
coefficients Cj and Cg.s- The original presentation of this phenomenological potential model in 
Ref. |15j sought to provide a comparison with the lattice calculations of the time, which had 
an incomplete classification of the HLHL states in terms of the total isospin I and spin S of the 
system, while also maintaining a connection with the physical B meson states. Because of this, 
classification of the various potentials was made in terms of the physical B and B* (first angular 
excitation of the B meson) with respect to the quantum numbers / and S. 

The difference in angular momentum spaces of the non-static and static limit prevents a direct 
interpretation of the lattice data from the present work in terms of physical B and B* states, and 
our classification of states makes it difficult to reconcile the previous classification with ours. We 
therefore choose to recalculate the spin dependent coefficients of the potential model relevant for 
the static limit BB system we study on the lattice, the results of which are presented in Table [n| 
(For details of the calculation, see Appendix [C]) . The previous determination of these coefficients 
for the HLHL system included spin degrees of freedom for the heavy quarks in the two meson states 
\MiMj > allowing for better classification of the potential in terms of non-static limit states. We 
choose to neglect the spin degrees of freedom of the heavy quarks in our determination, effectively 
fully implementing the static limit for the the potential model. Thus the spin degrees of freedom 
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(l,l^,|JJ,P^ PP.) = (0,0,0,+,-,+) Fit 




0.8 1 
d (fm) 



FIG. 5: Fit of the potential model in eq. [20|to the {I,I^,\J^\,Pj_,P,Pi) = (0,0,0,+,-,+) channel. The 
colored band represents the uncertainty in the fit paraniters (3 and g from jackknife analysis. 

of our two meson kets \MiMj > are just those of the spin of the hght degrees of freedom of our 
HLHL state. The evaluation of these coefficients however requires knowledge of the total angular 
momentum of the two meson state, a point that has been neglected until now. As we seek to fit 
the {I, Iz,\Jz\, P±, P, Pi) = (0,0,0,+,—,+) and (1,1,0,—,—,+) states, we need to determine if 
these particular states are in a symmetric angular momentum triplet, or an antisymmetric angular 
momentum singlet. In order to make this identification, we must rely on the overall symmetries of 
the state in question. We know that the parity P of a given state is the product of the intrinsic parity 
Pi and the symmetry of the spatial wavefunction. From this relationship, and with knowledge of the 
symmetry of the isospin spatial wavefunction, we can infer the symmetry of the angular momentum 
wavefunction: 



Symj = (-) iSymi) (Pi) (P) 



(22) 



where Symj and Symj the symmetries of the angular momentum and isospin wavefunctions. The 
overall negative sign appears from exchanging fermions in the parity operation. Using this we are 
able to identify the {I, Iz,\Jz\, P±, P, Pi) = (0,0,0,+,—,+) channel with Symj = — as a J = 
state, and the (/, Iz, \ Jz\, P±, P, Pi) = (1, 1, 0, — , — , +) channel with Symj = + as a J = 1 state. 
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The spin dependent coefficients can then be recalculated for our states and are shown in Table nU. 



C. Fitting Procedure and Bound State Determination 



In fitting the potential model of eq. 20 to our lattice data, we use two free fit parameters: /3 and 
g and take the remaining parameters b, rh and as to be O.lSGeV^, 0.33GeV, and 0.5 respectively 
as in Ref. [15]. A fit is performed for each of 305 jackknife ensembles, allowing for an accurate way 
to estimate the error on the extracted fit parameters, shown in Table [llj. As we are ultimately 
interested in the energy levels allowed by the potential model, and not the model parameters 
themselves, we will only briefly comment on the fit parameters. It is immediately obvious that 
g is not well determined for the J = 1 channel. It's also interesting that the fit parameter /3 is 
significantly smaller for the J = channel, indicating a much narrower spatial distribution of the 
two meson wavefunction. 

Once the fit parameters have been extracted they are then inserted into the two body radial 
Schrodinger equation to determine if any bound states exist. As we are restricting ourselves to 
L = states, the two body Schrodinger equation can be written as: 



2m dr^ 



+ F^"^' (r) 



u (r) = Eu (r) (23) 



where m is the reduced mass of a two B meson system (with the single meson mass taken from 
the Particle Data Group [16j), u (r) = (r) and V^^^ (r) is the potential model presented in the 
preceding section excluding the image terms. 



Eq. 23 is then solved numerically as an eigenvalue problem with a spatial discretization of 0.01 



fm and a spatial cutoff of 10 fm (corresponding to a sphere with r = 10 fm), and the boundary 
condition that ^ (r) = 0. This spatial volume provides ample space for the potential to decay 

r=10 

to zero. The eigenvalue spectrum is then analyzed for each of the two states discussed above. 
While the J = 1 channel exhibits a near continuum of positive eigenvalues (discrete only because 
of the numerical solution method), the J = channel does admit a single bound state with energy 
Eq = —50.0(5.1) MeV (with the uncertainty determined by carrying through the jackknife analysis 



from the fit parameters and solving eq. 23 for each of the 305 {(3,g) sets). Aside from the binding 
energy, we can also calculate the RMS radius for the two meson wavefunction ^ (r) from the 
wavefunctions u (r) above: 
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2\i/2_ E^rt\n{nr 

2 



(24) 



For the bound state wavefunction uq (r), we find an RMS radius of 0.383(6) fm, the error again 
estimated by jackknife analysis. 

Although no previous calculation of the binding energy in this particular static-limit channel 
exists (lattice or otherwise), Ref. |3| does quote binding energies and RMS radii for a doubly bot- 
tom (L, 5, /) = 0^ (0, 0, 0) channel which is consistent (in the static limit) with the quantum 
numbers of our static limit (/, Iz, | J2I, -P, Pj) = (0, 0, 0, -|-, — , -|-) channel. This reference uses 
two different potential models to calculate binding energies: the constituent quark cluster model 
CQC and the the Bhaduri-Cohler-Nogami or BCN model. The BCN model includes the same 
interactions as those used in Ref. |15j to derive the potential used to fit our lattice results (namely, 
color coulomb, linear confinement and spin-spin). Furthermore, the BCN parameters corresponding 
to string tension 6, strong coulpling as, and constituent quark mass fh used in [4j are very sim- 
ilar to those used in our potential model (compare our (6, as,m) = (0.18 GeV^, 0.5, 0.33 GeV) 
to (0.186 GeV^ ,0.52,0.337 GeV) ) . These binding energies should provide a relevent point of 
comparison for our results provided our lattice discretization errors have minor effects on the 
extracted potential model fit parameters. In comparison, we find our values for the binding 
energy and RMS radius to be consistent with the values quoted in [4J from the BCN model 
{Eq^thms) = (— 52MeV, 0.334fm), providing a good cross check that our lattice calculation has 
identified a bound state in the static limit (/, J^, | J^l, Pj_, P, Pj) = (0, 0, 0, -|-, — , -|-) channel. The 
fact that the bound state identified in that work has an RMS radius that is smaller than the 
sum of the individual mesonic RMS radii is indicative of the compact nature of that bound state. 
Additionally, as illustrated in Ref. [H] (see eqns. 4), the static limit HLHL tetraquark state can 
be written as a linear combination of products of two single meson wavefunctions in different spin 
states. This is consistent with the idea that although the compact tetraquark state may have a 
complicated color space structure composed of color vectors, this state can always be decomposed 
into a linear combination of products of two single meson wavefunctions. 



VII. CONCLUSIONS 



We have computed using lattice QCD the interaction potential between two b-meson states 
in the limit of static b quarks. With this lattice potential parametrized with a functional form 
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{J, J.) 


Ci 


Css 


/3 (GeV) 


9 


xVd.o.f. 


Eq (MeV) 


(0,0) 


-1 


3/4 


0.274(14) 


0.041(12) 


0.9943 


-50.0(5.1) 


(1,0) 


1 


1/4 


0.459(38) 


0.016(20) 


0.4119 


N/A 



TABLE II: Spin dependent coefRcients from reference [T^ and fit parameters from fitting our lattice data 
to a modified version of the model presented in ref. [151 . Here /3 corresponds to the spatial width of the HL 
meson wavefunction, and g is the coupling strength of the additional Yukawa term introduced in this work. 
The uncertainties quoted for the fit parameters are jackknife estimates. 

motivated by the quark model description of the two b-meson interaction, we have determined 
the bound state energies in the heavy-light-heavy-light (HLHL) tetraquark system. To perform 
this study we introduced colorwave propagators for calculating meson correlation functions and 
extended the formalism to the HLHL system in order to provide a novel way for an efficient 
calculation of HLHL correlation functions for several {I,Iz,\Jz\ , P±, P, Pi) channels. The effect 
of limiting the colorwave plane wave basis on the ground state overlap of single HL correlation 
functions was explored, and a choice for the momentum cutoff p^^^ was made to optimize the 
quality of the signal versus the computational cost. For a single HL meson, results indicate that a 
more localized interpolating field has a better overlap on the ground state, suggesting the compact 
nature of the HL meson. 

HLHL potentials were calculated for 24 distinct (/, J^, \ Jz\ , P±, P, Pi) channels, exhibiting three 
distinct asymptotic values as r — t- oo corresponding to the different ways B and Bi mesons can 
be combined. The tendency of the HLHL energy to overshoot the expected asymptotic value of 
+ Eb and 2Ebi may be due to contamination from excited states and the possibility of Bi 
mixing with a B — tt state. It was determined that the attractiveness or repulsiveness of the HLHL 
potential corresponds directly to the symmetry of the two meson spatial wavefunction under spatial 
inversion, in agreement with Ref. [H]. The asymptotic behavior of the various HLHL states was 
shown to be dependent on the intrinsic parity of the state. While the Pj = — states have only 
one asymptotic value (corresponding to a single two meson BBi component), the Pj = + channels 
have two asymptotic values corresponding to both BB and BiBi two meson components. By 
examining the construction of single HL correlation functions, it was determined that we could 
increase overlap with the BB and BiBi two meson wavefunctions by projecting the correlation 
functions to include only positive or negative parity components of the Dirac basis quark spinors. 

The existence of bound states was then explored for the (I, I^, \ Jz\, P±, P, Pi) = (0, 0, 0, +, — , +) 
channel as it exhibited a wider and deeper potential when compared with the other attractive 
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potentials. Analysis was also carried out for the {I, Iz,\Jz\, P±, Pi) = (15 1,0,—,—,+) for the 
purposes of comparison. A modified version of the potential model described in Ref. [15j was used 
to fit the lattice data, and two fit parameters /3 (the gaussian width of the HL meson wavefunction) 
and g (the Yukawa interaction constant) were extracted from, the fit. Inserting the potential with 
the extracted fit parameters into the two body Schrodinger equation, we then solved numerically for 
the eigenvalues of the hamiltonian, searching for any negative energy eigenstates. A single negative 
energy bound state was found in the (0, 0, 0, +, — , +) channel, with an energy of Eq = —50.0(5.1) 
MeV and RMS radius trms = 0.383(6) fm. These results were found to be consistent with results 
presented in Ref. [3j for the state (L, S, I) = 0+ (0, 0, 0) (which maps onto our (0, 0, 0, +, — , +) 
channel in the static limit). The errors quoted on these results are statistical only. One needs 
to account for several systematic errors such as 1/mf, corrections (m^ the b quark mass), lattice 
spacing effects as well as dependence on the light quark mass. ^ 
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Here we show that correlation functions for our B{Bi) states are composed entirely of up- 
per(lower) components of the Dirac basis components of the light quark flavors. We begin with a 
general HL correlation function with arbitrary source and sink operators (neglecting color indices 
and working in the Dirac basis): 
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Appendix A: Parity content of HL interpolating operators 



X 



= Y,{Q (x, t) T,q (x , t) q (x, 0) F^Q {x, 0)> 



X 




(Al) 



X 



Just before completion of this manuscript a study of the same system appeared as a preprint 18 . Some of these 
systematics were studied there. 
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Sh is a heavy quark propagator given by: 



Sh {x, t- to) = i \ W (f , t- to) = P+W {x, t; to) 




(A2) 



where W {x, t; to) is a Wilson Hne from to to t. Substituting this, we have: 



Chl (t) 



tr (75 (p+^t (^^ t- 0)) ^^t.Sl (x, t; 0) T,) 



2? 



^ tr, (W^ (x, t; 0) tra {TiP^TjSL {x, t; 0))) 



(A3) 



X 



where we have used 75P+75 = P_. For Fj = Tj = 1, we project to only the lower components 
of the Dirac basis light quark propagator, while for Fj = Tj = 75 we project only to the upper 
components of the Dirac propagator. 



To determine two quark wavefunctions in spin and flavor space yielding the quantum numbers 
(I, Iz,\Jz\, P±, Pi Pi), we begin with states of definite (I, Iz, J, Jzi Pi)'- 



where m,t,p are the projections of spin and isospin along the z-axis and the intrinsic parities 
of the light quarks, and the wi'i'^ma = (1/2, mi, 1/2, 7^2! J, J^), = {l/2,ti,l/2,t2\I, h) are 

the Clebsch-Gordon for angular momentum and isospin. From these operators, we average over 

= ±1 states and determine P± from the quantum numbers Pi and P and the spatial symmetry 
of the operator. It should be noted here that there are two combinations of (^1,^2) that contribute 
to the Pi = +1 HLHL states, and we make the decision to keep these as distinct operators. 

Linear combinations of the above operators are then taken to produce states of definite exchange 
parity P, the necessary combinations determined by summing over sets of the above operators that 
map onto each other under P with the appropriate weight Wp^ p^ = ±1 



Appendix B: Construction of light quark wavefunctions 





mi ,77*2 
tl,t2 



(Bl) 



(i,h,\j,\,p^,p,p,) 



Yl ^Lp2 tei (1^1)92(^2)] 



(B2) 



Appendix C: Determination of spin coefficients for potential model 



Here we present our derivation of the spin coefficients (7/ and C^s presented in Table [11} In Ref. 
[T7] . an interaction potential for two meson states is calculated by including spin-spin, color coulomb 
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and linear confinment interactions in a two quark interaction hamiltonian. By considering these 
interactions between each of the quark quark pairs in a 4 quark (2 meson) scattering state, transfer 
matrix elements are calculated and then Fourier transformed to give a corresponding position space 
potential. In Ref. [15j, this method was applied to the HLHL system. When calculating the spin 
dependent portion of the potential, all but one of the interaction diagrams (referred to as "Transfer 
2" ) can be neglected because the spin of the heavy quarks is neglected in the static approximation. 
This diagram includes an insertion of the interaction hamiltonian between the two light quarks, 
as illustrated in in Fig. [6| The spin dependent contribution of this diagram to the potential can 
be factorized such that all the dependence enters through two coefficients, which are defined as 
matrix elements between the initail and final two meson states: 



Ci = {CD\I\AB) (CI) 
Cs.s = {aD,\Si-Sj\AB,) (C2) 

Where / here is understood to be the identity operator in spinor space. Upon inspection of the 
diagram, it's clear that the matrix element of / will not always trivially be unity due to the quark 
interchange between the initial and final two meson state. 

With respect to Fig. |6l these matrix elements as outlined in [T7] are defined explicitly as: 



Cj = {CD\I\AB) 



Sb,si 



Cs-s = {CiDj\Si ■ Sj\AiBj 



(C3) 



(C4) 



For our purposes, we wish to entirely neglect the spin of the heavy quarks in the above matrix 
elements. Because of this, the Clebsch-Gordan coefficients XscSc (relating the spin of the two 
quark state to the meson state) are all unity. The states between which we wish to calculate these 



matrix elements are two particle angular momentum eigenstates | J, Jz)a h= ^b) , of which 
we are only interested in |1,0) and |0,0). To account the light quark exchange in Fig. [6| we note 
the following relations: 



\Sa,Sb) 



J=0, 

Jz=0 



1 

71 



(in) - lit)) 



1 



V2 



(lit) - Iti)) 



\Sc,Sd) 



.7=0, 

Jz=0 



(C5) 
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B 



D 



and 



\Sa,Sb) 



J=l, 

Jz=0 



FIG. 6: Transfer 2 diagram from Ref. 



^ (Iti) + lit)) = ^ (lit) + Iti)) = \Sc, S,) 



J=l, 

Jz=0 



(C6) 



From the above relations, it is easy to calculate the matrix elements of interest for our problem 
(for the states |1,0) |1,0) and |0,0) |0,0)): 



(i,o|,,,i|i,o),,, = (i,oL_,i|i,o)„,, = i 

(0, 0|,,, 1 10, 0),,, = (-) (0, ou,, 1 10, 0),,, = -1 



and 



(C7) 
(C8) 



(1, 0|,,, S^ ■ Sj |1, 0),_, = (1, 01^, S. • S, |1, 0),,, = 1/4 

(0, 0|, , S, • S, |0, 0) . , = - (0, 0|„ , S, • S, |0, 0)„ , = -(-3/4) 



(C9) 
(CIO) 
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